CERN-PH-TH/2006- 198 



QCD with light Wilson quarks on fine lattices (I): 
first experiences and physics results 

L. Del Debbio 1 , L. Giusti 2 '*, M. Liischer 2 , R. Petronzio 3 , N. Tantalo 3 ' 4 

1 SUPA, School of Physics, University of Edinburgh, Edinburgh EH9 3JZ, UK 

2 CERN, Physics Department, TH Division, CH-1211 Geneva 23, Switzerland 

3 Universitd di Roma "Tor Vergata" and INFN sezione "Tor Vergata", 
Via della Ricerca Scientifica 1, 1-00133 Rome, Italy 

4 Centro Enrico Fermi, Via Panisperna 89 A, 1-00184 Rome, Italy 
Abstract 

Recent conceptual, algorithmic and technical advances allow numerical simulations of lat- 
tice QCD with Wilson quarks to be performed at significantly smaller quark masses than 
was possible before. Here we report on simulations of two-flavour QCD at sea-quark masses 
from slightly above to approximately 1/4 of the strange-quark mass, on lattices with up to 
64 x 32 3 points and spacings from 0.05 to 0.08 fm. Physical sea-quark effects are clearly seen 
on these lattices, while the lattice effects appear to be quite small, even without 0(a) im- 
provement. A striking result is that the dependence of the pion mass on the sea-quark mass 
is accurately described by leading-order chiral perturbation theory up to meson masses of 
about 500 MeV. 



1. Introduction 



Many different formulations of lattice QCD are currently in use. The aim to reduce 
the lattice effects and to preserve chiral symmetry as much as possible has been the 
principal motivation for the introduction of increasingly complicated lattice actions. 
Highly improved actions are not obviously the best choice in practice, however, since 
they tend to slow down the numerical simulations by a large factor. Moreover, the 
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conceptual transparency of the lattice theory (otherwise one of its greatest assets) 
may be compromised in extreme cases. 

The philosophy advocated here is to keep the theory simple at the fundamental 
level (thus sticking to Wilson's formulation [1] and its close relatives) and to develop 
adapted computational strategies that allow simulations of large lattices with small 
spacings to be performed efficiently. Extrapolations to the continuum limit will then 
still be required, but the hope is that, in many cases of interest, the lattice effects 
will already be small at the accessible lattice spacings. 

For many years, the Wilson theory had the reputation of being difficult to simulate 
at light-quark masses significantly smaller than half the strange-quark mass. The 
situation has now changed completely, following the development of the DD-HMC 
simulation algorithm [2-4] and of a fine-tuned version of the Hasenbusch-accelerated 
HMC algorithm [5-8], both being much faster than the algorithms used thus far (for 
related earlier studies of unquenched lattice QCD, see refs. [9-14], for example). 
The success of these algorithms is partly also due to the fact, discovered later [15], 
that the Wilson-Dirac operator has a safe spectral gap in the large-volume regime 
of QCD, even though chiral symmetry is not exact in the Wilson theory. 

In this paper we report on extensive simulations of two-flavour QCD in a range 
of lattice spacings, lattice sizes and quark masses not explored before. Apart from 
sect. 2, where the simulations are briefly described, the emphasis is put on the main 
physics results. Further details of the simulations, data tables, etc., will be published 
in a forthcoming more technical paper. 



2. Simulation table 

As already indicated, the lattice theory is set up following Wilson [1], with a doublet 
of sea quarks of equal mass. In this theory, the sea quarks represent the up and down 
quarks, while the strange quark will be added later, at the level of a valence quark, 
i.e. without the associated quark determinant. 

Three series of lattices were simulated, labelled A, B and D (see table 1) f. In the 
D series, the Sheikholeslami-Wohlert term [16] was included in the quark action with 
coefficient c sw set to the value determined by the ALPHA collaboration [18], thus 
ensuring non-perturbative on-shell 0(a) improvement. The size of the representative 

| Our notation and normalization conventions coincide with those of ref. [17]. As usual we quote 
the values of (3 = 6 /</q and n = 8 + 2mo instead of the bare coupling go and sea-quark mass mo • 
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Table 1. Simulations included in the physics analysis 



Run 


Lattice 


P 






K 


iVcfg 


A, 


32 x 24 3 


5.6 





0, 


15750 


64 


A 2 








0. 


15800 


109 


A 3 








0. 


15825 


100 


B 1 


64 x 32 3 


5.8 





0. 


15410 


100 


B 2 








0, 


15440 


101 


B 3 








0. 


15455 


104 


B 4 








0. 


15462 


102 


D x 


48 x 24 3 


5.3 


1.90952 


0. 


13550 


104 


D 2 








0. 


13590 


171 


D 3 








0. 


13610 


168 


D 4 








0. 


13620 


168 


D 5 








0. 


13625 


169 



ensemble of gauge-field configurations generated in each case is given in the last 
column of table 1. 

In physical units, the lattice spacing on the A, B and D lattices is estimated to be 
0.0717(15), 0.0521(7) and 0.0784(10) fm respectively (see sect. 4). The three series 
of simulations cover a similar range of sea-quark masses, from values slightly above 
the strange-quark mass m s down to values close to m s /4. 

All simulations were performed using the DD-HMC algorithm [4] , which combines 
domain decomposition ideas with the Hybrid-Monte-Carlo algorithm [19] (hence the 
name DD-HMC). As explained in refs. [4,21], the speed of the algorithm also very 
much depends on the use of the Sexton- Weingarten multiple-time integration scheme 
[20]. In practice, a relevant performance figure is the number 7V op of floating-point 
operations required for the generation of an ensemble of 100 statistically independent 
field configurations on a 2L x L 3 lattice at a specified lattice spacing and sea-quark 
mass. A formula that fits our experience with the DD-HMC algorithm well is 

,/20MeV\/ L V/CUfmV^n . , 

where m denotes the running sea-quark mass in the MS scheme at renormalization 
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Fig. 1. Quark-line diagrams contributing to the two-point correlation function (3.1). 
In the second diagram, two pairs of sea-quarks are created from the vacuum and bind 
into a pair of pseudo-scalar mesons. 

scale /i = 2 GeV and k ~ 0.05 if the 0(a)-improved theory is simulated (k ~ 0.03 
without improvement). 

In 2001, at the annual conference on lattice field theory in Berlin, a similar formula 
was presented by Ukawa [22], summarizing the experience made by the CP-PACS 
and JLQCD collaborations with the algorithms available at the time. With respect 
to that formula, the scaling exponent of the quark mass in eq. (2.1) is reduced from 
3 to 1, the exponent of the lattice spacing from 7 to 6, and the constant k is roughly 
100 times smaller. 

Apart from the operations count, the suitability of the simulation algorithm for 
parallel processing is a key issue, if large lattices are to be simulated. Domain decom- 
position methods tend to perform well from this point of view, and one of the design 
goals of the DD-HMC algorithm was in fact to keep the communication overhead 
small [3,4]. Special-purpose computers are then not required and most simulations 
(including the B series) were actually performed on commodity PC clusters with up 
to 64 double-processor nodes f. 



3. Physical sea-quark effects 

Once the sea quarks are included in the simulations, an obvious question is whether 
their presence has a visible effect on the computed quantities. Correlation functions 
of local fields depend on the sea-quark content of the theory in various ways. The 
flavoured pseudo-scalar densities, for example, can only couple to multi-meson states 
through the creation of virtual quark pairs (see fig. 1). Higher-states contributions 
to the correlation functions of these fields are therefore expected to increase when 
the sea-quark mass is lowered (thus moving away from quenched QCD). 

| The program that was used for the simulations of the O(a)-improved theory (the D series) can 
be downloaded from http://cern.ch/luscher/DD-HMC. 
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Fig. 2. Simulation results (data points) for the effective pseudo-scalar mass M e f[(t) 
in MeV as a function of the time t in fm (runs D2 and D4). The average valence-quark 
masses m va i and the meson masses M (grey bands) are nearly the same in the two 
cases, while the sea-quark mass m sea changes by a factor of 2 (quoted quark masses 
are bare current-quark masses). 

Indications for the presence of multi-meson intermediate states in some two-point 
correlation functions were already found some time ago by the UKQCD collaboration 
[23]. We now have data at smaller quark masses, where the effect is more pronounced 
and where a significant dependence on the sea-quark mass is seen in both the pseudo- 
scalar and vector channels. 

For illustration we introduce two valence quarks, labelled r and s, and consider 
the two-point correlation function 



C(t) 



d 3 x ((r 75 s) (x) (s 75 r) (0)) 



(3.1) 



x =t 



of the corresponding flavoured pseudo-scalar density at zero spatial momentum (for 
simplicity we use a continuum notation in this section). In finite volume, C(t) may 
be expanded in a spectral series 



C(t) = c e- Mt + Cl e- M,t + ... 

t — >oo 



(3.2) 



where M denotes the mass of the associated pseudo-scalar meson and M' the energy 
of a three-meson state with all particles at rest. 
Plots of the effective mass 



M eS (t) = -— lnC(t) = M + ce" 



-(M'-M)t 



+ ... 



(3.3) 



now show that the higher-states contributions are not small in general. Moreover, 
as is evident from fig. 2, they tend to grow when the sea quarks become lighter. In 
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the cases shown in the figure, the energy M' of the lowest three-meson intermediate 
state is expected to be approximately equal to M+2M W , where M n denotes the mass 
of the pseudo-scalar mesons made of the sea quarks. The two-state formula (3.3) 
actually fits the data quite well if this expression is assumed and if M n is determined 
from the sea-quark pseudo-scalar correlation function (solid lines in fig. 2). 

While the observed enhancement of higher-states contributions is in line with qual- 
itative theoretical expectations, their presence also tends to complicate the analysis 
of the simulation data. In particular, at small sea-quark masses, the computation 
of hadron masses may require accurate data at larger time separations than was the 
case in quenched QCD. Multi-mass fits and variational methods can be helpful at 
this point, although the associated systematic uncertainties must then be balanced 
against the possibly lower statistical errors. 



4. Setting the scale 

The choice of a physical reference scale is an important step in the analysis of the 
simulation data. Results obtained on different lattices can then be expressed in units 
of this scale and thus be compared with one another. 

Following ref. [17], we adopt a mass-independent scheme where the lattice spacing 
in physical units is the same on all lattices at a given bare coupling. Different choices 
of the reference scale are possible, none of which appears to be free of some practical 
or conceptual shortcoming. Here we add a valence strange quark to the theory and 
determine the scale through the pion mass and the masses of the pseudo-scalar and 
vector mesons that are made of a strange antiquark and a sea quark (we refer to 
these as the K and the K*). More precisely, we adjust the quark masses so that the 
ratios Mk/Mk* and M^/Mk assume some prescribed values and then take Mk as 
the reference scale. 

Since we wish to set the scale in a physically sensible way, we require the ratio 
Mk/Mk* to be equal to its physical value of 0.554. This condition fixes the strange- 
quark mass m s at any given coupling and sea-quark mass m (see fig. 3). Ideally we 
would like the latter to be such that M^/Mk assumes its physical value too, but 
this would require a long extrapolation in the sea-quark mass and, moreover, would 
be a point where the K* is unstable (i.e. the extrapolation would have to go through 
a kinematical threshold). 

We now note, however, that once m s is fixed, the reference scale in lattice units, 
ciMk, appears to be weakly dependent on m, particularly so at small m (see fig. 3). 
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Fig. 3. On each lattice, the masses of the K and the K* were computed at 4 or 
5 values of the bare strange-quark mass. The results obtained on the D lattices are 
plotted here for illustration (data points). At the lighter sea-quark masses, the point 
where Mk/Mk* = 0.554 (dotted line) is found by a quadratic interpolation in the 
strange-quark mass of the nearest data points. The solid line shows the interpolation 
in the case of the lattice D$. 



The reason for this behaviour (which is seen on all series of lattices) could be that 
both Mk and Mk* are functions of m + m s rather than of m and m s separately, up 
to corrections proportional to the squares of the masses. In any case, the observation 
suggests the reference scale to be defined at the point where, say, M^/Mk = 0.85, 
which is within the available data range. This convention, although somewhat un- 
physical, is entirely satisfactory for the purpose of comparing results from different 
lattices. 

The results for the reference scale obtained in this way are summarized in table 2. 
In order to avoid any confusions, we mark all quantities evaluated at the reference 
point with a subscript "ref " . The sea-quark and strange-quark hopping parameters 
at the reference point, for example, are denoted by K re f and K s , re f- Setting Mx, re i = 
495 MeV, this leads to the lattice spacings quoted in the last column of the table, 
while for the pion masses at the smallest sea-quark masses on the A, B and D series 
of lattices we obtain 403, 381 and 377 MeV respectively. 

The lattice spacings calculated here are significantly smaller than those previously 
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Table 2. Determination of the lattice spacing* 



Lattice series 


K ref 


^s,ref 




a [fm] 


A 


0.15822(3) 


0.15769(4) 


0.180(4) 


0.0717(15) 


B 


0.154561(12) 


0.154257(10) 


0.1310(17) 


0.0521(7) 


D 


0.136207(7) 


0.135912(13) 


0.197(3) 


0.0784(10) 



* At the quark masses where M K /M K * = 0.554 and M n /M K = 0.85 



published by us in a conference report [21], where the Sommer radius [24] was used 
as reference scale. Larger lattice spacings are also obtained if the scale is set by 
the K and K* masses, similarly to what was done here, but at larger sea-quark 
masses (see fig. 3). As a result of the new determination of the lattice spacings, our 
estimates of the pion masses in MeV are pushed to higher values than those quoted 
in ref. [21]. Moreover, we decided to discard the simulation at the lightest quark 
mass reported there, because the lattice turned out to be too small for that mass. 

All this illustrates the fact that at present the assignment of physical units remains 
somewhat ambiguous. For a definitive solution of the problem, simulations at smaller 
quark masses will probably be required, and the scale setting may eventually have 
to be based on the properties of the stable hadrons (the pion and the nucleon in the 
two-flavour theory). 



5. Quark-mass dependence of M n and F n 

The bare quark masses that appear in the lattice action of the Wilson theory require 
a power-divergent additive renormalization. This complication can be bypassed by 
extracting the quark masses from the PCAC relation (see ref. [17], for example; in 
the improved theory, we used the non-perturbatively improved axial current [25]). 
Moreover, ratios of these masses do not need to be renormalized since the multi- 
plicative renormalization factor cancels. 

In fig. 4 the ratio (M^/M^ref ) 2 is plotted as a function of the corresponding ratio 
of quark masses. If there were no systematic effects, all data points shown in this 
figure would have to lie on a single curve, within statistical errors, representing the 
mass dependence of {M n /Mx^ei ) 2 in the continuum limit. Note that the statistical 
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Fig. 4. Dependence of the square of the pion mass M T on the sea-quark mass m. 
The solid curve is a quadratic least-squares fit (with constant term) of all data points, 
and the plot on the right is a blowup of the region enclosed by the little box. 



errors of the points are uncorrelated, except for the correlations that are introduced 
through the normalization factors. The quality of the empirical fit (solid line) then 
suggests that no statistically significant systematic effects are seen in this plot. 

Another rather striking outcome is that M% is a nearly linear function of the sea- 
quark mass m in the range covered by the data. There is a visible curvature towards 
the larger masses in fig. 4, but the coefficient of the quadratic term in the empirical 
fit, y = -0.03(3) + 1.03(5)x + 0.02(2)x 2 , is small. In the range M n /M K , ief < 1.1, 
the data are also well represented by a straight line through the origin. 

The corresponding plot of the pion decay constant F x , given in units of the decay 
constant Fx,tei of the K meson at the reference point, is more difficult to interpret 
(see fig. 5). Apart from the fact that the statistical errors tend to be larger here, the 
results of the D series of simulations appear to be significantly different from those 
of the A and B series. There is no obvious curvature in either set of data points, 
and correlated straight-line fits (solid lines) are found to be statistically consistent. 
Although the two lines are visibly different, the fitted values of their slopes, 0.235(11) 
and 0.192(11), deviate from each other by less than 3 times the combined statistical 
error. The statistical significance of the effect is thus not overwhelming. 

When discussing systematic errors, an important point to note is that the axial- 
current renormalization constant cancels in the ratio of decay constants plotted 
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Fig. 5. Dependence of the pion decay constant on the sea-quark mass m. The 
solid curves are linear least-squares fits of the data points from the A and B lattices 
(upper line) and of the points from the D lattices (lower line). 



in fig. 5. Moreover, the ratio is largely insensitive to the values of the improvement 
coefficients ca and 6 a [17,25] on which the axial current in the 0(a) improved theory 
depends. Lattice effects may still be significant, however, and we can also not exclude 
the presence of important finite-volume effects. A more specific problem is that any 
variations in the normalization factors Fx^ei an d m re f + Tn s>xe i propagate to the 
slopes of the lines in fig. 5. The alignment of the data points may therefore appear 
to be better or worse, depending on the statistical fluctuations at the reference point 
and on its detailed specification. 

The fact that the lines in fig. 5 have nearly the same intercept in the chiral limit is 
probably an accident. Both lines also practically pass through F^/Fx^ef = 0.82 (the 
experimental value of F^/Fx) at M^/Mx ref = 0.28. However, as will be shown in 
the next section, such extrapolations to smaller quark masses could be misleading. 



6. Comparison with chiral perturbation theory 



In two-flavour QCD with unbroken isospin symmetry, the chiral expansion of the 
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pion mass reads [26] 



M 4 

Ml = M 2 + MM 2 /A|) + . . . , M 2 = 2Bm, (6.1) 

where F, B and A3 are a priori unknown constants. A phenomenological analysis, 
taking low-energy experimental data as input, suggests [26,27] 

F = 86.2 ± 0.5MeV, l 3 = ln(A 2 /M 2 )\ M=139 6MeV = 2.9 ± 2.4, (6.2) 

while B depends on the renormalization scheme for the quark mass and thus cannot 
be determined from such data alone. The chiral expansion of the pion decay constant 
has the form [26] 

M 2 

\n{M 2 /A 2 ) + ... (6.3) 



l6n 2 F 

and the phenomenological discussion leads to the estimate [28] 

h = HA 2 /M 2 )\ M=139 6Mey = 4.4 ± 0.2 (6.4) 

for the low-energy constant A 4 . 

In principle the low-energy constants can be determined from lattice data without 
recourse to phenomenological estimates. However, as will become clear shortly, the 
available simulation data are insufficient for a solid analysis of this kind. While the 
situation in the case of the pion mass is somewhat more favourable, our principal 
goal in the following will be to find out whether the data are compatible with the 
expansions (6.1) and (6.3) for a reasonable choice of the parameters. 

We first need to rewrite the equations in a form where all dimensioned quantities 
are expressed in units of the scales at the reference point. To this end, it is helpful 
to introduce the abbreviations 

m ref + m s>ref ' 327r2i? i, re f ' 

and to define the scaled parameters 

F = ~E?—i B= m ^ 2 ms ^ B, l n = HKl/M 2 KM ). (6.6) 

^K,ref M K,ref 
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2m/(m ref +m jref ) 



Fig. 6. Fit of the quark- mass dependence of the square of the pion mass M, in 
the range M^/Mx^ef < 1-1, using the one-loop formula (6.7) with C — 0.072 and 
F — 0.70 (solid line). The fit is a correlated least-squares fit of all data points in the 
specified range, with unconstrained parameters B and I3. 



The chiral expansions 



M 2 ~ B 2 r 2 
_E_ = Bx + C {\ n (Bx) -l 3 } + ..., (6.7) 

M K,ref F z 



Br 

2C^-{\n(Bx)-h } + -.., (6i 



FR,ref F 

may now be directly compared with the simulation data (note that l n = l n — 2.53 if 
M K ,r e { = 495 MeV is assumed) . 

The computation of the decay constant Fx^ef, an d thus of the constant C, in- 
volves the renormalization constant Za of the axial current. Recent estimates of the 
latter in the two-flavour Wilson theory at the couplings of the A and B lattices are 
0.77(2) and 0.78(2) [29], while in the case of the D series of lattices we may use the 
value 0.75(1) determined by the ALPHA collaboration [30]. For the constant C we 
then find 0.068(4), 0.071(4) and 0.076(3) respectively. These figures are barely con- 
sistent with one another, suggesting the presence of lattice or finite-volume effects, 
but it should also be noted that the determination of Za is not free of systematic 
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Fig. 7. Fit of the quark-mass dependence of the pion decay constant F v in the range 
M n /M K ,ref < 1-1, using the one-loop formula (6.8) with C = 0.072 and B = 1.106 
(solid line). The dashed line with its 1-sigma error margin (grey band) represents an 
alternative fit that includes a hypothetical two-loop term. 

ambiguities [29,30] f. 

A very accurate determination of C is fortunately not needed for the chiral fits, 
because C only appears at next-to-leading order in the chiral expansions. We thus set 
C = 0.072 and simplify the fit procedure by substituting F = 0.70 in eq. (6.7), which 
will turn out to be an approximately correct value. In the range M^/Mx^ef < 1.1, 
the one-loop formula (6.7) then fits the data for the pion mass very well, the fit 
parameters being B = 1.11(6)(3) and I3 = 0.5(5)(1) (see fig. 6; the second errors are 
estimates of the systematic uncertainty arising from the inaccurately known values 
of C and F). We did not attempt to estimate the effects of any higher-order terms in 
eq. (6.7) so that the quoted values of the fit parameters should be taken as effective 
values, describing the data in the specified range of pion masses. 

In the case of the pion decay constant, the comparison of the simulation data with 
the chiral formula (6.8) is complicated by the scattering of the data points in fig. 5, 
which may partly be the result of systematic effects. However, since the points line 

I The same comments apply in the case of Fx^ei where we obtain 107(3), 105(3) and 101(2) MeV 
on the A, B and D lattices (assuming M^^ef = 495 MeV as before). The fact that these results are 
all lower than the decay constant Fx = 113(1) MeV of the physical kaon may not be significant, 
because the two-flavour theory neglects the effects of the strange sea quark. 
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up at the smaller quark masses, we may attempt to fit these, setting C and B to 
the previously determined values and adjusting F and I4 (see fig. 7). The statistical 
quality of this fit (solid line) turns out to be quite good, but the curvature of the 
fit function is not seen in the data and the fit therefore appears to be somewhat 
artificial. 

A more plausible fit (dashed line) can be obtained by including a hypothetical two- 
loop term proportional to B 2 x 2 /F 3 in the chiral expansion (6.8), with a coefficient 
C = 0.046 that is not unreasonably large. The fit parameters F and Z4 change from 
0.60(4) and 1.6(1) to 0.73(3) and 0.73(8), respectively, when the two-loop term is 
added. 

The discussion in this section shows that simulation data at significantly smaller 
quark masses, with small systematic and statistical errors, will be required for a 
reliable determination of the parameters in the chiral lagrangian. It seems safe to 
conclude, however, that our results in the range M n /Mk, ref < 1-1 are not incom- 
patible with chiral perturbation theory. In particular, the fact that M 2 is a nearly 
linear function of the quark mass m is not in conflict with the presence of the chiral 
logarithm in eq. (6.7). 



7. Concluding remarks 

In the coming years, simulations of lattice QCD with Wilson quarks will no doubt 
rapidly progress towards smaller quark masses and lattice spacings than are reported 
here. In order to guarantee the stability of the simulations [15], but also to keep the 
finite-volume effects under control, the constraints 

M„L>3, L>2fm, (7.1) 

should be respected in these computations. On a given lattice, the bounds (7.1) set 
a lower limit on the lattice spacings and pion masses that can be reached (see fig. 8). 
Simulations of the two-flavour theory may not be practical at all these points, but 
the cost formula (2.1) is encouraging and suggests that simulations at a < 0.08 fm 
and M n < 300 MeV, for example, can be performed already with the computer 
resources available at present. However, to be able to sort out the systematic errors, 
many lattices will have to be simulated which may require a coordinated effort. 

On the A, B and D series of lattices, the smallest values of M^L are 3.5, 3.2 and 
3.6 respectively, while the spatial sizes L of the lattices are estimated to be 1.72, 1.67 
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Fig. 8. Range of lattice spacings a and pion masses M w denned by the bounds (7.1) 
on a 2L x L 3 lattice as a function of the lattice size L/a (shaded area above the line 
labelled by the corresponding value of L/a). 

and 1.88 fm, i.e. somewhat below the required minimum. Finite-volume effects may 
not be totally negligible on these lattices and will need to be investigated, extending 
the studies by Orth et al. [14] to smaller quark masses and lattice spacings. So far 
we did not include the nucleons in the physics analysis, because these are probably 
even more sensitive to finite- volume effects than the mesons. 

It may be somewhat surprising that no significant lattice effects were seen in fig. 4, 
even though O(a) counterterms were only included in the D series of simulations. 
The weak dependence on the lattice spacing could be related to the fact, first noted 
by Sharpe and Singleton [31], that the O(a) lattice effects amount to an additive 
quark-mass renormalization at leading order of chiral perturbation theory. Since the 
quark masses that appear in the PC AC relation already include all additive renor- 
malizations, it follows that the data points plotted in fig. 4 are insensitive to these 
leading-order lattice effects. 

The mass dependence of the pseudo-scalar decay constant, on the other hand, is a 
second-order effect in chiral perturbation theory. At this order, only some of the O(a) 
lattice corrections can be compensated by a renormalization of the parameters in 
the chiral lagrangian, and an accidental O(a) improvement is therefore not expected 
in this case. 

We wish to thank Rainer Sommer for technical discussions and Gilberto Colangelo, 
Stephan Diirr, Jiirg Gasser and Heiri Leutwyler for some very helpful notes, summa- 
rizing some relevant results of chiral perturbation theory. The numerical simulations 
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